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We study the finite temperature properties of two-component fermionic atoms trapped in a two- 
dimensional optical lattice. We apply the self-energy functional approach to the two-dimensional 
Hubbard model with a harmonic trapping potential, and systematically investigate the thermody- 
namic properties of this system. We find that entropy and grand potential provide evidence of a 
crossover between the Mott insulating and metallic phases at certain temperatures. In addition, we 
find that entropy exhibits a cusp-like anomaly at lower temperatures, suggesting a second or higher 
order antiferromagnetic transition. We estimate the antiferromagnetic transition temperatures, and 
clarify how the trapping potential affects this magnetic transition. 

PACS numbers: 03.75.Ss, 05.30.Fk, 67.85.Lm, 75.30.Kz 
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I. INTRODUCTION 

Ultracold atoms in an optical lattice created by inter- 
fering laser beams are bridging the gap between theoret- 
ical and experimental studies in the field of condensed 
matter physics. We can control the depth, interstice and 
dimension of the lattice by manipulating lasers, and the 
interaction between ultracold atoms by using the Fesh- 
bach resonance [1|. This high degree of controllability 
allows us to simulate correlated many-body systems ex- 
perimentally 0t3 • Theoretically, it has been pointed out 
that this system can be regarded as a realization of the 
Hubbard model, which is one of the most fundamental 
models including the many-body effects [a-S]- I^ f^^t, 
various phenomena described by this model have been 
successfully demonstrated in experiments |9l4l5l|. 

The first experimental progress was reported as an 
observation of the quantum phase transition of bosonic 
^^Rb atoms between a superfluid and a Mott insulator in 
an optical lattice [9]. As regards the fermionic atoms, the 
Fermi surface and its topological change were observed 
in the fermionic ^°K optical lattice system, suggesting 
a metal to band- insulator transition [lOJ. Recently, by 
combining numerical and experimental studies, it was 
clarified that a Mott transition was realized in ''"K opti- 
cal lattice systems |li|,[i3]- In these reports, the Hubbard 
model with a trapping potential was analyzed by using 
a local density approximation with the dynamical mean- 
field approach (LDA-I-DMFT) [l^] and a zero-tunneling 
limit calculation [Ilj. The effects of the trapping po- 
tential, which is an important characteristic of ultracold 
atoms in an optical lattice, have been properly taken into 
account in addition to those of the correlations. However, 
these calculations present certain difficulties in terms of 
investigating magnetic transitions. The observation of 
the magnetic ordered phase in the optical lattice system 
is a major concern for condensed matter physicists, be- 
cause it could provide ways to elucidate the nature of high 
Tc superconductors [l^ . By using the recently developed 
real-space dynamical mean-field theory (R-DMFT) [iTl - 



|20|, it has been pointed out that an antiferromagnetic 
(AF) ordered phase is stable in the Hubbard model with 
a trapping potential at zero temperature 18]. This nat- 
urally motivates us to undertake a detailed study of both 
the magnetic transition and the Mott transition at finite 
temperatures using a reliable numerical method. 

For this purpose, we investigate the two-component 
Fermi-Hubbard model in a two-dimensional (2D) opti- 
cal lattice with a harmonic trapping potential at finite 
temperatures. We make use of the self-energy functional 
approach (SFA), which has been successfull y appli ed to 
homogeneous Hubbard- type lattice models |2l| - [24| . We 
demonstrate that this method can properly take account 
of the effects of the trapping potential in addition to those 
of strong correlations and thermal fiuctuations. SFA fur- 
ther provides us with important thermodynamic quanti- 
ties such as grand potential and entropy. We find that 
both Mott and AF transitions can be characterized by 
the behavior of these thermodynamic quantities. The 
AF transition temperatures of the present model are sys- 
tematically examined by varying the temperature, inter- 
action strengths and curvatures of the harmonic trapping 
potential. 

This paper is organized as follows. In Sec. [Hi we in- 
troduce the two-component Fermi-Hubbard model on a 
2D optical lattice with a harmonic trapping potential. In 
Sec. IHIl we briefiy outline the application of SFA to the 
present system. In Sec. IIVI we discuss the finite temper- 
ature properties of the Mott insulating region. In Sec. |Vl 
we discuss how the AF ordered region is realized in the 
system, and estimate the AF transition temperatures. In 
Sec. IVIl we briefiy summarize this paper. 



II. MODEL 

In this paper, we investigate the two-component Hub- 
bard model on a 2D lattice with a trapping potential 
that we assume to be harmonic. The Hamiltonian of this 
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FIG. 1: (Color online) A schematic diagram of the 2D lattice, 
where the cross is the center of the harmonic trapping poten- 
tial, L is the total number of sites, a is the lattice distance, 
and ri is the distance between the center of the trap and the 
jth site. The numbers written in the top right of the lattice 
sites correspond to the site indices. 



sulating phases |T21, is defined as, 

Rsc = {r)/rt, 



0- wE^'^^/^^^' 



(4) 

(5) 



where iV^ — '^^{nia) is the number of atoms at the ith 
site, and A'tot — X]j ^i i^ the total number of atoms. 
To discuss the magnetic ordered phase, we calculate the 
magnetization Mi = {ni^^) — {ni{). We also investigate 
thermodynamic quantities: grand potential 51 and en- 
tropy S — —dft/dT, which are sensitive to both the Mott 
transition and the AF transition. 



model is given by H == "Ht + "Hu, 

^t = -^ E E4^j- + EE(^*'-?-^)"-' (1) 

■Hu = y^^UnifUji, (2) 



where cl^icia) creates (annihilates) a fermionic atom 
with pseudospin (7{—^,^) at the ith site, and Uia ~ 
cl^Cia- Here, we describe the nearest- neighbor hopping 
integral as J, the curvature of a harmonic trapping poten- 
tial as Vt, the chemical potential as /^, and the interaction 
strength between two atoms with different pseudospins as 
U. For simplicity, we define the site index i as shown in 
Fig. [TJ where L is the number of sites, a is the lattice 
distance, and r^ is the distance between the ith site and 
the center of the trap. It is useful to introduce the pa- 
rameter matrix t, which characterizes the non-interacting 
Hamiltonian "Ht as 

Vtrf iii ^ j 

[t]^^. = \ J if the ith site neighbors the jth site (.3) 
otherwise 



We define the scales of length and energy in the follow- 
ing way such that the calculated results do not depend on 
the details of the parameters, e.g. L,a,Vt, and fi. We in- 
troduce the characteristic trap length rt = \/Na-/y/Tra, 
which corresponds to the radius of a non-interacting 
atomic cloud in the zero-tunneling limit, where N^ is the 
total number of atoms with pseudospin a. The distance 
ri is rescaled by this n as rsc = ri/rt- We also introduce 
the bandwidth W = 8J as a. scale of energy. We consider 
the characteristic trap energy Et = Vtr^, which corre- 
sponds to the Fermi energy in the zero-tunneling limit 
p^ . In addition to Et/W, we consider two parameters: 
the temperature T/W and the interaction strength U/W. 

Before concluding this section, we summarize the phys- 
ical quantities calculated in this paper. The rescaled 
cloud size, which characterizes the Mott and band in- 



III. METHOD 

The thermodynamic properties of the present sys- 
tem are studied with the self-energy functional approach 
(SFA) [HHll], which is based on the Luttinger-Ward 
variational method [2^. It has been pointed out that 
SFA allows us to undertake an efficient investigation 
of the finite temperature properties of homogeneous 
Hubbard-type lattice systems, for instance, the infinite- 
dimensional Hubbard model [2ll . |23 | and also the 2D 
Hubbard model [2J|, taking account of the effects of 
strong correlations. We extend this method to deal with 
inhomogeneous systems. Here, we explain an applica- 
tion of this SFA to the Hubbard model with a trapping 
potential. 

First, we begin by briefly explaining the general frame- 
work of the SFA. According to the Luttinger-Ward func- 
tional approach [25], the grand potential fi of a Hubbard- 
type system described by the Hamiltonian H = Ht + Hu 
is written as, 



r2[S] = F[S] + Trln[(Go- 



(6) 



where F[E] is the Legendre transformation of the 
Luttinger-Ward potential, S is the self-energy, and Go = 
(w — /x — 1)~^ is the non-interacting Green function. Here, 
we use the notation TrA = T^^ ■[A]--{iuJn), where 
ujn = (2n + l)TrT is the Matsubara frequency. Un- 
der the condition 9ri[S]/9S = 0, we obtain the phys- 
ical Green function G that satisfies the Dyson equation 
G^^ = Gq"^ — S. In general, the functional F[S] is not 
known explicitly, which prevents an evaluation of 57 [S]. 
However, the potential i^[S] does not depend on the de- 
tails of the non-interacting Hamiltonian "Ht as long as 
the shape of the interaction term 'Hu remains unchanged 
[2l| . This allows us to introduce a reference system that 
has a Hamiltonian with the same interaction term as that 
of the original system. This reference Hamiltonian is ex- 
plicitly given by 'H'^°^ — Tit' + T~i\] with the variational 
parameter matrix t'. The grand potential of the reference 
system 17'''^^ is then written as, 
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F[S]+Trln[(Go 



'-1 



(7) 



where Gq = (w — /x — t')~^. By subtracting Eq. ([7]) 
from Eq. ^ , we can rewrite the grand potential of the 
original system as a function of the self-energy for the 
reference system S''°^: 

n[S''°f] = f7''"f + Trln[(w + M-t-S''°f)-i] 

- Trln[(w + /i-t'-S'-'=f)-i] . (8) 

A reference system with an optimized parameter matrix 
t' satisfying the condition 



dn[Y,"'^]/dt' = 0, 



(9) 



gives us an appropriate self-energy S''°^, Green function 
G = {uj-fi-t- i:'"^)-^ and grand potential r2[S'°f], 
which approximately describe the physical quantities of 
the original system. 

An application of the SEA to the Hubbard model with 
a trapping potential is achieved as follows. The simplest 
reference system with which to investigate the present 
model is L-sets of two-site Anderson impurity models 
[22| . Here, the ith impurity site, corresponding to the 
ith site of the original lattice, is connected to the ith 
non-interacting atomic bath. The Hamiltonian of this 
reference system is given by the form: "H™* = X]i^i°^j 
and 

a 

+ Y^ {y^„4^a^„ + H.c^ + Uui^n,^, (10) 

(T 

where ci[„{airT) creates (annihilates) an atom with pseu- 
dospin a at the ith atomic bath. The variational pa- 
rameter Via is the hybridization of the impurity and the 
atomic bath, e^ is the effective potential, and hi is the ef- 
fective magnetic field. Here, we briefly explain the role of 
variational parameters e^, hi and Via- By optimizing an 
effective potential ti under the condition dfl/dej = 0, the 
number of atoms iV^ is properly adjusted [23|. We can 
discuss the stability of magnetic ordered phases via the 
condition dQ/dh, = ^^,^. The hybridization V,a 
effectively describes the hopping integral between the zth 
site and adjacent sites in the original lattice; therefore, 
we can also discuss the Mott transition via the conditions 
dn/dV,a = [H [13. Details of the role of the atomic 
bath with the hybridization term will be discussed in ap- 
pendix A. 

By employing exact diagonalization, we can easily ob- 
tain the grand potential ilp^ and the self-energy S^J^^ of 
the ?th reference system. Now, the grand potential of the 
original system is given by. 
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'TrlnGl'^f 
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Grcf //^rcf — 1 v'l'cf \ — 1 



(11) 
(12) 

(13) 



where [GJ^^JJ^j = Sij/{iuj, 
rewritten as. 



£» - ahi - Vil/iujn), and 
The variational condition ^ is 



TV([G.],,-[Gf], 



9s: 



ref 



dt' 



0, 



(14) 



where we denote the variational parameter sets as t^ = 
{e,,h,,V,^,Vi^}. Note that [G-^]^,^ and d^f/dt[ de- 
pend only on one parameter set t'j. In contrast, [Ga-]jj 
depends on all of the elements of the variational param- 
eter matrix t', and we have to solve an equation of the 
4Lth degree. It is still difhcult for the large size sys- 
tem to solve equation p4p and to optimize the varia- 
tional parameter matrix t'. To avoid this difficulty, we 
introduce another parameter matrix t* and local Green 



functions as GfJ'{t',t*) 



(g: 



cav— 1 / 



^ZHt')): 



where Gr(t*) = (s™f(t*) + [G.(t*)]^::^j . Were- 

place [Gct]^^ in equation (fT4|) with G\'^ considering the 
additional self-consistent condition [Ga]-- — G^°^. If t* is 
given, the variational condition (|14p is finally rewritten 
as. 




T^lnG^; 



T^lnGl- =0. 



(15) 



We note that the equation of the 4ith degree is decom- 
posed into independent L-sets of the 4th degree. Now, we 
choose t* as a parameter matrix satisfying this new con- 
dition ()15p , and repeatedly solve this variational problem 
until t* converges (at the same time, the self-consistent 
condition [Gct],;j = G'°'^ is satisfied). In summary, we 
self-consistently solve independent L-sets of easily solv- 
able SEA problems instead of a huge problem with the 
4Lth degree. 

The procedure mentioned above can be extended to 
any other inhomogeneous systems. In addition, if we 
deal with the attractively interacting systems, we can 
straightforwardly extend the scope of this method to the 
s-wave superfluid phase by adding the term, AicLcL -I- 
H.c, to the reference Hamiltonian (ITU]), where A^ is a 
variational parameter corresponding to the superfluid or- 
der parameter. Eurthermore, the SEA has been suc- 
cessfully applied to multi-component Eermi systems [28J , 
where it has been suggested that novel quantum phase 
transitions could be observed |29l - [33 |. Promising candi- 
dates for such systems are fermionic ^^^Yb atoms |14| 
and/or ^Li atoms [Ullil. 

Now, we describe briefly how to calculate the physical 
quantities. It is possible to evaluate various quantities 
denoted in the previous section by means of the Green 
function G(aj), the grand potential fi, and its deriva- 
tives. Eor instance, the total number of atoms is ob- 
tained by the relations: A'tot = dVl/dfi or A^tot — TrG. 
Note that the variational condition dft/dt' — guaran- 
tees that these two relations are equivalent to each other 
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FIG. 2: (Color online) Rescaled cloud size Rsc as a function 
of Et/W for different U/W at T/W = 0.1. 



[35| . In addition, the SFA allows us to calculate quan- 
tities in momentum space, e.g. experimentally observed 
time of flight images. We can also calculate angle re- 
solved photo emission spectra using the present scheme, 
which is closely related to a recent observation by the 
JILA group in Ref. [36|. However, it is beyond our cur- 
rent scope to investigate these quantities in momentum 
space. 

We close this section with few comments on the supe- 
riority of the SFA to other related numerical approaches. 
The zero-tunneling limit calculation and LDA+DMFT 
present certain difficulties when we discuss quantum 
phase transitions accompanied by spontaneous symmetry 
breaking. Although R-DMFT can be applied to the rnag- 
netic ordered phase 18] and the superfluid phase J19l.l20|. 
it is difficult to investigate the finite temperature prop- 
erties systematically and precisely with this method. In 
contrast, the SFA provides comprehensive ways to inves- 
tigate various quantum phase transitions at finite tem- 
peratures. Moreover, this method allows us to calculate 
various useful quantities observed experimentally. 



IV. MOTT INSULATING PHASE 

In this section, we investigate the finite temperature 
properties of the 2D Hubbard model with a harmonic 
trapping potential. In particular, we focus on the be- 
havior of the thermodynamic quantities of this model at 
certain temperatures, where the Mott insulating region 
is observed. 

First, to discuss the validity of our method, we cal- 
culate the rescaled cloud size Rsc which was previously 
evaluated for the three dimensional (3D) model by the 
LDA-hDMFT approach in Ref. ^. Figure H shows 
Rsc as a function of Et/W for different U/W values, 
where the temperature is fixed at T/W — 0.1. When 
the interaction strength is small {U/W = 0, 0.5 and 
1.0), Rsc decreases rapidly in the weakly trapped re- 




2 

1.5 

fe 1 

0.5 




(c) E/W~l 


(d) E/W~2 


^^^ 


"'"'^'^ 


^x/Vv 


- ^\ - 




. ^k-^ 



0.5 



1 
''sc 



1.5 



0.5 



1 

''sc 



1.5 



FIG. 3: (Color online) Density profiles of atomic cloud, i.e. 
the number of atoms Ni as a function of the rescaled distance 
rsc for several fixed values of Et/W: (a) 0.6, (b) 0.8, (c) 1.0 
and (d) 2.0. The interaction strength is varied from U/W = 
0.5 to U/W = 2.0. The temperature is fixed at T/W = 0.1. 



gion Et/W < 1, while it changes little in the strongly 
trapped region Et/W > 1. It should be noted that, in 
the strongly trapped limit Et/W :$> 1, Rsc finally con- 
verges to l/-\/2 ~ 0.7, which corresponds to the cloud size 
of non-interacting atoms without tunneling. This indi- 
cates that most of the atoms except for those around the 
site ri ~ rt become inactive; in other words, the band- 
insulating like state appears [l2] • On the other hand, for 
large U/W values of 1.5 and 2, the Rsc curves exhibit 
shoulder like structures in the region 0.4 < Et/W < 0.7 
and 0.4 < Et/W < 1.0, respectively, suggesting the for- 
mation of the Mott insulating state there. 

To clarify these points further, we calculate the density 
profiles, namely the number of atoms Ni as a function of 
rescaled distance rsc = fi/ft, for several fixed values of 
Et/W and U/W. We show the results in Fig. [3] Since we 
adopt a grand canonical ensemble in our calculations, the 
Et/W values in Fig. [3] have a small deviation of about 
±0.05. In Fig. |3i;a), at Et/W - 0.6, a weU developed 
Mott plateau with Ni = 1 from the center to the edge of 
the cloud appears for large U/W values of 1.5 and 2.0. 
These plateau profiles are gradually deformed as Et/W 
is increased. From Fig. ^h) and (c), we see that the 
additional atoms are loaded in the Mott plateau around 
the center at Et/W ~ 0.8 {Et/W ^ 1.0) when the inter- 
action strength is U/W = 1.5 (2.0). In Fig. Efc) and (d), 
the Mott plateau is not dominant and is limited to the 
shell region around rsc ~ 1 even for the large interaction 
U/W = 2.0. In aU the panels in Fig. O for smaU U/W 
values of 0.5 and 1.0, there are no Mott plateaus over 
the entire Et/W range. This density profile behavior is 
consistent with the shoulder structure of the Rsc curves 
in Fig. El 

Our results for the 2D system exhibit good qualitative 
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FIG. 4: (Color online) (a) Entropy per atom S/Ntot and (b) 
grand potential per atom Q/Ntot for several values of U/W 
at T/W — 0.1. Thin (light blue) lines in the panel (b) mean 
Q — —EtNtot/2 and Q = —EtNtot from top to bottom. 



agreement with previous theoretical and experimental re- 
sults for the 3D system reported in Ref. [l^]- This sug- 
gests that a proper renormalization allovirs us to compare 
the results between two- and three-dimensional optical 
lattice systems qualitatively. 

As a complementary study, we next investigate 
other thermodynamic quantities at a fixed temperature 
T/W = 0.1. We calculate the entropy per atom S/Ntot 
and the grand potential per atom Q/Ntot as shown in 
Fig. mja) and (b), respectively. Curves of S/Ntot reveal 
a similar Et/W dependence of Rsc to that in Fig. [2] For 
strongly interacting cases such as U/W =1.5 and 2.0, the 
shoulder structure appears in the S/Ntot curves around 
Et/W ~ 0.5, suggesting that the strong correlations in- 
duce localized free spins whose entropy takes a constant 
value of In 2. It is reasonable to expect that a magnetic 
ordered phase will appear because these free spins are 
interacting with each other. However, at T/W = 0.1, 
thermal fluctuations will destroy such an ordered phase, 
as we discuss in the next section. 

Grand potential also provides us with useful informa- 
tion on the Mott transition of the present system. As 
shown in Fig. IDJb), when U/W — 1.5 (2.0), the gradient 
of Q/Ntot curves gradually changes around Et/W ~ 0.75 
(1.0) indicating that there is the crossover between the 
Mott and metallic phases. Here, the Mott (metallic) 
phase is defined as the phase in which the Mott insu- 
lating (metallic) region is dominant. We find that, in the 
Mott and metallic phases, the grand potential obeys the 
relations fi ex —EtNtot and ex — i?t7Vtot/2, respectively. 
Note that the characteristic trap energy Et corresponds 
to the energy required to add one atom around the edge 
of the atomic cloud, namely the chemical potential in 
the limit of zero-tunneling. The difference between these 
relations, a factor 2, results from the fact that double 
energy is required to add one atom in the Mott phase. 

We now shift our attention to the effects of thermal 
fluctuations on the quantities shown above. In Fig. [SJ we 
show the temperature dependence of the Rsc curves for 
U/W = 1.5. Figure [S] shows the temperature dependence 
of the density profiles for U/W — 1.5 a.t Et/W ^ 0.6 and 
1.0. 




FIG. 5: (Color online) Rescaled cloud size Rsc vs Et/W for 
different T/W values. The interaction strength is fixed at 
U/W = 1.5. 



1.5 



(a) E/W~Q.6 



(b) £/W~0.8 



r/W=0.30 
0.20 X 
0.15 K 
1-1 »»-«,. 0-10 

•I, 0.05 • 



0.5 



iiM 






0.5 1 1.5 0.5 1 1.5 2 

''SC ''SC 

FIG. 6: (Color online) Density profiles of atomic clouds for 
two characteristic trap energies: (a) Et/W ~ 0.6 and (b) 
Et/W ~ 0.8. The temperature is varied from T/W = 0.05 to 
T/W = 0.3. The interaction strength is fixed at U/W = 1.5. 



From Fig. El we find that, below T/W = 0.10, Rsc 
curves have a shoulder structure around Et/W ~ 0.5, 
and change little as T/W is decreased. On the other 
hand, above T/W = 0.10, thermal fluctuations destroy 
the shoulder structure; therefore, we find no signature 
of the Mott insulating state in the Rsc curves. How- 
ever, as shown in Fig. [HI the Mott plateau survives 
up to T/W ^ 0.15, suggesting that the destruction 
of the shoulder structure of Rsc curves is mainly at- 
tributed to the thermal excitation of metallic atoms at 
the edge of the Mott insulating region. Consequently, for 
Et/W ^ 0.6 and 0.8, we roughly estimate a specific tem- 
perature T* /W value of ~ 0.15 around which the Mott 
plateau is destroyed by thermal fluctuations. 

We also calculate the entropy per atom S/Ntot by vary- 
ing the temperature as shown in Fig. [71 The S/Ntot 
curves do not saturate at low temperatures {T/W < 0.1) 
in contrast to those of Rsc- At the lowest temperature 
T/W = 0.025 in Fig. [71 the entropy per atom is highly 
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FIG. 9: (Color online) Real-space distribution of the magne- 
tization at U/W = 1.5 and T/W = 0.025 for {a.)Et/W ~ 0.6 
and {h)Et/W ~ 0.8. 



FIG. 7: (Color online) Entropy per atom S/Ntot vs Et/W 
for several T/W values. The interaction strength is fixed at 
U/W = 1.0. 
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FIG. 8: (Color online) (top) the magnetization Mi and (bot- 
tom) the number of atoms Ni as functions of the rescaled dis- 
tance rsc for two characteristic trap energies: {a.)Et/W ^ 0.6 
and {h)Et/W ~ 0.8. The interaction strength is varied from 
U/W — 0.5 to U/W = 2.0 and the temperature is fixed at 
T/W = 0.025. 



suppressed in the region 0.4 < Et/W < 0.8. This is 
due to the fact that a magnetic ordered phase appears at 
T/W = 0.025, as discussed in detail in the next section. 
We note that entropy is sensitive to magnetic transitions. 
As shown in this section, SFA allows us to calculate 
various thermodynamic quantities at finite temperatures. 
Such thermodynamic quantities have already been inves- 
tigated experimentally in ultracold atomic gases trapped 
in a conventional magnetic trap |37h39| . We hope that 
these thermodynamic quantities will be measured in op- 
tical lattice systems in the near future. 



V. MAGNETIC ORDERED PHASE 

In the previous section, we showed that entropy shows 
a characteristic of the AF magnetic transition at a low 
temperature T/W = 0.025. Here, to clarify the prop- 
erties of the magnetic ordered phase, we investigate the 
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FIG. 10: (Color online) Absolute value of the magnetiza- 
tion \Mi\ as a function of rsc and T/W for fixed values of 
{a)Et/W ~ 0.6 and {h)Et/W ~ 0.8. 



present model at a fixed temperature T/W — 0.025. 

We first look at the magnetization Mi = (n^^) — (n^). 
In Fig. Elja) and (b), we show Mi and Ni as functions of 
rsc — Ti/rt for different V /W values at Et/W ^ 0.6 and 
0.8, respectively. Except for a weakly interacting case of 
U/W = 0.5, we find finite values of M,. The sign of Mi 
changes alternately from site to site, suggesting the AF 
order. To clarify this point, we also show the real-space 
distribution of the magnetization Mi for U/W = 1.5 in 
Fig. [9] We note that the rotational symmetry is broken 
in the AF ordered region. The magnetization is well- 
developed in the Mott plateau region. In addition, even 
in the metallic region in the vicinity of the Mott plateau, 
the magnetization Mi has a small but finite value [Ig . 

From Figs. [5] and [51 we see that the density profiles of 
atoms remain unchanged at very low temperatures even 
though the magnetization Mi becomes finite. As a con- 
sequence, the appearance of a magnetic ordered region 
hardly affects the rescaled cloud size Rsc (see also Fig. 

El). 

The magnetic ordered phase has already been investi- 
gated at zero temperature using R-DMFT by M. Snoek 
et al. in Ref. [l8|. Our results at very low temperature 
agree well with their results at zero temperature. We 
next focus on the temperature dependence of magneti- 
zation, and evaluate the critical temperatures of the AF 
magnetic transition in the present system. 

We systematically calculate the magnetization by vary- 
ing the temperature. Figure [TUf a) and (b) show the ab- 
solute values of Mi for U/W = 1.5 at Et/W - 0.6 and 
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FIG. 11: (Color online) Entropy per atom S/Ntot as a func- 
tion of T/W for several values of Et/W at U/W — 1.5. Inset 
is an enlarged view and arrows indicate the magnetic transi- 
tion temperatures. 



0.8, respectively. The magnetic ordered region gradu- 
ally spreads as T/W decreases. Thus we can understand 
the quantum phase transitions in the present system as 
follows. The Mott plateau first develops at higher tem- 
peratures T*/W{^ 0.15) (see Fig. [5]), and then the mag- 
netic ordered region gradually grows from the inside of 
the Mott plateau as the temperature decreases below the 
transition temperature Tc/W{^ 0.04). Note that a well- 
developed Mott plateau region {Et/W ~ 0.6) has higher 
transition temperatures than those of the shell- like region 
(Et/W ^0.8). 

For a more precise discussion, we calculate the temper- 
ature dependence of the entropy per atom. In Fig. [TTl 
we show S/Ntot as a function of T/W at U/W — 1.5. 
For Et/W - 1.2,0.8,0.6 and 0.3, the S/Ntot gradient 
suddenly becomes steeper at the transition temperatures 
TJW = 0.029,0.038,0.048 and 0.031, respectively. In 
the region from just above the transition temperatures 
Tc/W to the specific temperatures T* /W, we find a linear 
temperature dependence of S/Ntot oc T. Significantly, 
except for the small Et/W value of ^ 0.3, the curves of 
the entropy per atom show a cusp-like anomaly at T^/W , 
suggesting second order magnetic transitions. However, 
we cannot definitely determine the order of the transition 
in our present calculation. We reach the conclusion that 
the transition is of the second or higher order. 

Let us discuss this temperature dependence of the en- 
tropy per atom in more detail. In a uniform system 
{Vt — 0), the linear temperature dependence of entropy 
5" ex T is a characteristic of a metallic state, while a con- 
stant entropy 5* = In 2 is a characteristic of the Mott 
insulating state at which a localized free spin is induced 
at each site. On the other hand, in a trapped system 



{Vt y^ 0), both states coexist below T*, namely the char- 
acteristic temperature of a well-developed Mott plateau; 
therefore, in the region Tc < T < T* , the entropy per 
atom obeys the relation S/Ntot — 6T-fc In 2, where b and 
c are, roughly speaking, the coefficient depending on the 
density profiles and the effective mass of atoms. As the 
temperature is decreased below T* , free spins in the Mott 
insulating region interact with each other, and then or- 
der antiferromagnetically at T = Tc where the redundant 
entropy cln2 is suddenly released. Above T*, a compli- 
cated temperature dependence can be seen that results 
from thermal fluctuations of the atoms in the Mott insu- 
lating region. Note that the entropy, which is enlarged 
by strong correlations, and its release play key roles for 
the magnetic transitions. Therefore, we first observe the 
developed Mott plateau around T*, and then find that 
the magnetic ordered region appears inside this plateau 
at Tc. 

Next, we systematically estimate the AF transition 
temperatures Tc/W for several choices of U/W = 1.0, 1.5 
and 2.0. The results are shown in Fig. [T^l In each 
curve, the transition temperatures as a function of Et/W 
have a maximum value around Et/W ^ 0.5 where we see 
the well-developed Mott plateau. When Et/W is small, 
there are insufficient atoms to form a magnetic ordered 
phase, leading to a very small Tc/W . The long tail of 
Tc/W curves in the large Et/W region results from the 
fact that the shell-like structure of the Mott plateau sur- 
vives in such regions. For comparison, we calculate T""': 
the AF transition temperatures for the uniform Hubbard 
model at half-filling {Vt — and iVtot = L)- We show 
T^'^-^/W for U/W = 1.0, 1.5 and 2.0 by the arrows in Fig. 
[T^ With a strongly interacting limit, T""' is inversely 
proportional to C/, while for a weakly interacting limit, 
T""' decreases exponentially with decreasing U [40| • In- 
deed, Tc^^/W has a maximum value around U/W ^ 1.0. 
On the other hand, in the system with a trapping po- 
tential, the transition temperatures for U/W = 1.0 are 
comparable to those for U/W = 1.5. This is because the 
Mott plateau does not develop at U/W = 1.0, because 
the effects of the confinement prevail against those of the 
correlations. For U/W = 1.5, we find higher transition 
temperatures than those for U/W — 1.0 and 2.0 over a 
wide Et/W range. The results in Fig. [T^suggest that the 
system for U/W ~ 1.5 is more suitable for an observation 
of the magnetic ordered phase. 

Before concluding this section, we compare our present 
results to other numerical works which have discussed 
the AF transition temperatures Tc in optical lattice sys- 
tems without consideration of the effects of the trap- 
ping potential [4l|, 14^. F. Werner et. al have estimated 
Tc/W ~ 0.04 at U/W ~ 1 [iil, and W. Hofstetter et. 
al have done Tc/W ~ 0.04 at U/W = 0.5 [4l|. Our es- 
timated value of Tc/W - 0.04 around Et/W - 0.5 at 
U/W = 1 shows a reasonable agreement with these pre- 
vious results. We note that, however, lower Tc values are 
realized except for around Et/W ^ 0.5 in Fig. [12] Fur- 
thermore, for a weakly interacting case U/W = 0.5, we 
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cannot find the AF transition down to T/W 
to the effects of the trapping potential. 
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FIG. 12: (Color online) AF transition temperatures Tc/W vs 
Et/W for different U/W values. The lines are guides for the 
eyes. The arrows indicate AF transition temperatures for the 
uniform Hubbard model at half-filling for U/W = 1.0, 1.5 and 
2.0 from top to bottom. 



d-wave superfluid correlations is not straightforward in 
our present choice of the reference system. These issues 
are beyond our current scope and left as important future 
works. 
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Appendix A: Role of the atomic bath 



(a) U/W=0.5 



(b) UAV=1.5 



1.5 

1 

0.5 



--. \^ 


V, -*- 


:::iH*^ -.- 






^=^ 


i. 


. \...... 



0.5 



1 

''sc 



1.5 



0.5 



VI. SUMMARY 

We have investigated the two-component fermionic 
atoms on a two-dimensional (2D) optical lattice with a 
harmonic trapping potential at finite temperatures. For a 
comprehensive understanding of both the Mott transition 
and the magnetic transition in this system, we have ex- 
tended the self-energy functional approach (SFA) to deal 
with inhomogeneous systems. By introducing additional 
self-consistent loops, a complicated variational problem 
for an inhomogeneous system following the framework 
of the SFA is decomposed into several easily solvable 
SFA problems. We have applied this method to the 2D 
Hubbard model with a trapping potential. A proper 
rescaling of the system parameters allows us to quali- 
tatively compare the two-dimensional system with the 
three-dimensional (3D) one. The calculated results of 
the rescaled cloud size and density profiles show good 
qualitative agreement with previous results in Ref. [12| . 
Furthermore, we have systematically calculated thermo- 
dynamic quantities such as entropy and grand potential. 
We clarified that entropy shows evidence of both the 
Mott transition and the antiferromagnetic (AF) transi- 
tion. In addition, we have demonstrated how confine- 
ment affects the AF transition temperature. We have 
estimated the AF transition temperature, and proposed 
a suitable parameter region for observing the AF ordered 
phase experimentally. 

A direct comparison with experimental results in the 
2D optical lattice systems will be useful. The detail infor- 
mation of experimental setups is required for such stud- 
ies. Additionally, it is very important to investigate d- 
wave superfluid phase |4l |. However, dealing with the 



FIG. 13: (Color online) The optimized hybridization Vi* as 
a function of rsc for two choices of the interaction strength: 
(a) U/W — 0.5 and (b) U/W = 1.5. Other parameters are 
fixed at T/W — 0.05 and Et/W ~ 0.8. For convenience, we 
also plot the number of atoms Ni. 

In this appendix, we comment on the role and impor- 
tance of the atomic bath with the hybridization Via in 
the reference Hamiltonian of Eq. ([TU]) . 

As mentioned in Sec. IIII[ we adopt the L-sets of two- 
site Anderson impurity models as the reference system. 
Here, we explain this procedure in more detail. We first 
separate a lattice into localized single sites. Then, for 
each separated site, we additionally introduce the site 
consisting of non-interacting atoms with the hybridiza- 
tion, and set the energy level of such sites just at Fermi 
level of the original system. Accordingly, these non- 
interacting sites are partially occupied and play a role 
of the atomic reservoir. Furthermore, the hybridization 
in the reference Hamiltonian corresponds to the hopping 
in the original Hamiltonian. 

It is known that, when the interaction strength is in- 
creasing, an effective energy scale of the hopping becomes 
much smaller. In other words, the renormalization effects 
are induced by the effects of correlations. We will demon- 
strate that the hybridization Via in the reference Hamil- 
tonian clearly refiects such many-body dynamics via the 
variational condition dft/dVia = 0. In Fig. [131 we plot 
the optimized values of the hybridization V*^ which sat- 
isfy all of the required conditions mentioned in Sec. IIIII 
Figure [iSlJa) and (b) show Vi* = V^ 



■n 



V*^ as a func- 



tion of rsc for U/W = 0.5 and 1.5, respectively. We 
set other parameters as T/W — 0.05 and Et/W '^ 0.8. 
In Fig. [T^Ta). we see large and constant Vi values in an 



occupied region {rsc ^ 1-3)- Note that, we reasonably 
choose 1^ = in an unoccupied region (rgc ^ 1-3) since 
the variational conditions dfl/dVia = is satisfied by 



any values of Via- On the other hand, in Fig. [T3l b). we 
see that values of Vi decrease in the Mott plateau region 
(0.5 < rsc < l)j indicating the renormalization effects. 
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